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Abstract 

We study the lattice reaction diffusion model 3A — > 4A, A — > ("triplet creation") using 
numerical simulations and n-site approximations. The simulation results suggest that the phase 
transition is discontinuous at high diffusion rates. In this regime the order parameter appears 
to be a discontinuous function of the creation rate; no evidence of a stable interface between 
active and absorbing phases is found. Based on an effective mapping to a modified compact 
directed percolation process, we shall nevertheless argue that the transition is continuous, despite 
the seemingly discontinuous phase transition suggested by studies of finite systems. 

PACS numbers: 05.70.Ln, 82.20.Wt 
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I. INTRODUCTION 



The exploration of phase transitions in simple, one component nonequilibrium models 
has attracted considerable interest 

Has, a, 

and important steps towards identifying 
the related universality classes have been achieved . In nonequilibrium models phase 
transitions may occur even in one-dimensional systems: the well known arguments, due 
to Landau and to van Hove |7|, against phase transitions in one-dimensional systems with 
short-range interactions, do not apply in the absence of detailed balance. However in low 
dimensions the effect of fluctuations is stronger, making continuous phase transitions more 
common. (A familiar example is the three-state Potts model, which exhibits a continuous 
transition in two dimensions, and a discontinuous one for d > 2.) 

In one dimension, discontinuous phase transitions have been found in models with long- 
range interactions j^l, or a conserved density 0, [h]], and in multi-component systems 
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141 ] . Compact directed percolation (CDP) has a discontinuous transition between a pair of 
absorbing states (all sites full or all empty) ; a similar transition between absorbing 

states is found in the one-dimensional Ziff-Gulari-Barshad model jl71 ]. A discontinuous 
transition between and active phase and an absorbing one in a single- component model was 
claimed for the triplet creation model (TCM) 18j, which does not possess a conservation 
law or long-range interactions. This model features particle reactions 3A — > AA, A — > and 
explicit diffusion (hopping) [19(. On increasing the diffusion probability, a crossover from a 
continuous to a discontinuous phase transition was detected in simulations and cluster mean- 
: ield approximations. Similar behavior was also reported in a stochastic cellular automaton 

Subsequently, Hinrichsen argued that in one dimension, discontinuous transitions between 



an active and an absorbing state cannot exist in models like the TCM 20|, |2l) . The original 



findings for the TCM [18( were nevertheless confirmed in spreading simulations by Cardoso 
and Fontanari 22] and in fixed order-parameter simulations by Fiore and de Oliveira 23]. 
The spreading exponents are shown in |22j to be those of compact directed percolation (CDP) 
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Park 



a trier itical point is suggested for a diffusion probability D ~ 0.95. Very recently 



241 ] reported simulation results that again support a continuous phase transition, 



belonging to the directed percolation (DP) universality class, at high diffusion rates. 

Recently, a field theoretic analysis of bosonic react ion- diffusion (RD) models led to a 
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hypothesis [5], based on a general phase transition classification scheme: bosonic, one- 
component RD systems with n-particle creation and m-particle annihilation always exhibit 

a first-order transition if n> m. This is indeed the case above the upper critical dimension 
I 

(see [25]). However in bosonic models one has to introduce a higher-order coagulation term 
m'A — > (m' — I) A with [m! > n), to avoid an infinite particle density in the active phase. 
Furthermore the topological phase space method used in deals with the reactions (creation 
and annihilation), but does not take into account the effect of diffusion, which turns out to 



be relevant in some cases, when different reactions compete |26|, |27|, 

In this work we study the TCM in an effort to determine whether multi-particle cre- 
ation, combined with rapid diffusion, can overcome fluctuations and generate a discontinuous 
phase transition in one dimension. This is a problem of longstanding interest in nonequilib- 
rium statistical physics, and is related to the existence of first-order depinning transition in 
nonequilibrium wetting (i.e., in a system with multiplicative noise, with an attractive wall) 
29]. 

The remainder of this paper is organized as follows. In Sec. II we define the model, and 
review applicable simulation methods and previous results regarding the nature of the phase 
transition. Section III is devoted to a discussion of n-site approximations, and Sec. IV to 
our simulation results. In Sec. V we use these results to motivate a simplified description of 
the model in the high diffusion rate regime, and discuss the nature of the transition using 
this mapping. Finally in Sec. VI we summarize our findings. 



II. THE TRIPLET CREATION MODEL 



The TCM is defined on a lattice, with each site either vacant or occupied by a particle; 
multiple occupancy is forbidden ■ In the one-dimensional TCM, a particle (A) attempts 
diffusion at rate D < 1, creation (3 A — > 4A) at rate A(l — D)/(l + A), and is annihilated 
(A — ► 0) at rate (1 — D)/(l + A). In a diffusion attempt, one of the nearest neighbor 
(NN) sites of the particle is chosen at random, and the particle jumps to this site if it is 
empty. If the target site is occupied, the configuration remains the same. In a creation 
attempt, if both NN sites of the particle are occupied, then one of the second- neighbor 
sites of the central particle is chosen at random, and if this site is empty, a new particle 
is placed there. If the conditions of two occupied NN sites and an empty target site are 
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not fulfilled, the configuration does not change. Annihilation occurs independently of the 
states of neighboring sites. The configuration with all sites empty is absorbing. Since the 
sum of these transition rates is unity, the total transition rate in a system with N particles 
is simply N. In simulations, the time increment associated with each attempted event 
(whether accepted or not) is At = 1/N, and one Monte Carlo step (MCS) corresponds to 
an accumulated time increment of unity. 



In Ref . 



18j the one-dimensional TCM was shown to exhibit a phase transition between 



the active and absorbing states; the transition was found to be continuous (and in the DP 
universality class) for smaller diffusion rates, but discontinuous for large D. By a discontin- 
uous transition we mean one in which the order parameter is a discontinuous function of the 
relevant control parameter(s), in the infinite-size limit. In the TCM the order parameter is 
the particle density p, and the control parameters are A and D. Since one of the phases is 
absorbing, at a discontinuous transition p should jump between zero and a finite value. 

The characterization of a transition as continuous or discontinuous in numerical simula- 
tions is fraught with difficulties: finite-size rounding can mask the discontinuity, and any 
finite system must eventually become trapped in the absorbing state. To circumvent these 
problems, a number of strategies have been proposed. 

Hysteresis with a weak source. A characteristic feature of discontinuous phase transitions 
is hysteresis. If one of the phases is absorbing, however, hysteresis cannot be observed simply 
by varying a control parameter, since the absorbing phase allows no escape. Bideaux, Boc- 
cara and Chate [30] showed that when the transition is discontinuous, adding a weak source 
of activity changes the absorbing and active phases to low-activity and high-activity phases, 
respectively. One may then observe a hysteresis loop between these phases, on varying the 



control parameter. This approach was used in [30] to demonstrate a discontinuous phase 
transition in a probabilistic cellular automaton, and was applied to the TCM in 18) , yielding 
a hysteresis loop. Below, we shall revisit the question of scaling under a weak source. 

Conserved order parameter simulations. In conserved order parameter simulations 311 ] . 
particles are neither created nor destroyed. Changes in configuration occur through particle 
jumps, which can be of any size up to that of the entire system, in a manner that respects 
the local rules of the process. The simulation yields an estimate for the control parameter 
value corresponding to the chosen order parameter density. Using this method, Fiore and 
de Oliveira found evidence for a discontinuous transition in both the TCM and the related 
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pair creation model (with creation reaction 2 A — > 3 A) at high diffusion rates 23]. 

Quasistationary (QS) simulation. As in conserved order parameter simulations, QS sim- 
ulation removes the absorbing state from the dynamics, but in a manner that samples the 
quasistationary probability distribution (i.e., conditioned on survival) 32]. A study of the 
TCM using this method 14| showed that as the system size tends to infinity, the QS order 
parameter appears to develop a discontinuity between zero and a positive value, as A is 
varied at a high diffusion rate, D = 0.98. (For a finite system the discontinuity is of course 
rounded.) A study of the TCM with biased diffusion (hopping in one direction only) yielded 



evidence of a sharp discontinuity 



33]. 



Spreading simulations. Studies of the spread of activity, starting from a seed at the origin, 
lave long been employed to characterize continuous phase transitions to an absorbing state 
l|, [34]. At the critical point, the survival probability P(t), mean number of active sites n(t), 
and mean-square distance R 2 {t) of active sites from the origin, all follow power laws. At a 
discontinuous transition, there is in principle no reason to expect scale-invariant spreading 
dynamics. Nevertheless, in the case of the TCM, Cardoso and Fontanari 22| demonstrated 
power-law spreading at the transition point, A c , for D = 0.98. The scaling exponents were 
identified as those of CDP, which, as noted above, suffers a discontinuous transition between 
a pair of symmetric absorbing states. 

Interface motion. Suppose we prepare the system with all sites occupied, allow it to relax 
to the QS state, and then remove all particles from half of the lattice. In the subsequent 
evolution, the interface between active and inactive regions broadens due to diffusion, and 
in general drifts toward one region or the other. Below we report studies showing that 
the drift velocity is proportional to A — A c . In a related analysis, we initialize the system 
with all sites in the region 1, ...,M occupied, and sites M + 1, ..,L empty, and study the 
long-time survival probability P(M). At the transition, the dependence of P(M) on M is 
consistent with independent, randomly diffusing interfaces, as in CDP. This result supports 
the existence of two phases, one absorbing, the other active, separated by a large gap in 
density. The two phases do not coexist: the fluctuating interfaces eventually meet, and one 
of the phases is lost from the system. 

Summarizing, the above mentioned studies, some from the recent literature, others to 
be reported below, provide evidence for a discontinuity in the QS order parameter, for 
hysteresis, and for a connection between the TCM at high diffusion rate and compact directed 
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percolation. 



A. Hinrichsen's objection 

Some years ago, Hinrichsen presented an argument to the effect that discontinuous phase 
transitions between an active and an absorbing state are impossible in one-dimensional 
systems with local interactions, and without additional conservation laws, special boundary 



conditions, or macroscopic currents 20, [2jJ. The argument is based on the observation 
that the effective surface tension of interfaces in such systems does not depend on the size 
of the domains they delimit. Hinrichsen's argument prohibits the presence of fixed, stable 
boundaries between coexisting phases; as noted, no such boundaries have been observed 
in simulations. But this in itself does not appear to imply that the dependence of the 
order parameter on growth rate must be continuous at the transition. The one-dimensional 
totally asymmetric exclusion process (TASEP), for example, exhibits a discontinuous phase 
transition in a certain region of parameter space, even though the position of the boundary 



between high- and .ow-density phases fluctuates over the entire system M . 

Hinrichsen |20J also reported simulation results supporting a continuous transition in the 
TCM at diffusion rate D = 0.9, that is, above the estimate for D t given in [l8[]. It is now 
generally acknowledged that D t > 0.9 in the TCM. 



In a recent study 24j , Park reported simulation results that support DP-like scaling in 
the TCM at diffusion rates 0.95 and 0.98. Specifically, the order parameter (starting from a 
filled lattice) appears to decay at long times as p(t) ~ t~ s , with S taking its DP value, over 
about two decades in time. We note however that the decay exponent is very sensitive to 
the choice of the time interval used for analysis and of the control parameter A. Analyzing 
simulation results for D = 0.98 in studies extending to 10 9 MCS, we obtain local decay 
exponents 8 e ff between 0.1 and 0.2, varying A in a very narrow range. A crossover between 
a long supercritical plateau for t < 10 7 and a rapid decay to an inactive state cannot be 



ruled out. For D = 0.98, the regime during which DP-like scaling is found in Ref. 24] 
(i.e., 10 7 <t< 10 9 ) corresponds to overall particle densities in the range 0.67 - 0.32. While 
scaling behavior can be observed at such densities in the contact process [l| , definitive results 
for the decay exponent would require studying systems with substantially smaller values of 
p. Finally, three or more critical exponents would have to be determined to demonstrate 
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convincingly that the transition falls in the DP class. 

Although the results of Hinrichsen and of Park do not appear to rule out rigorously a 
discontinuous transition in the TCM, we believe that they are fundamentally correct. This 
conclusion is based not on simulation results but rather on a mapping to a modified CDP 
process, to be developed in Sec. V. Analysis of this mapping in leads to the conclusion that 
the transition is in fact continuous, despite abundant numerical evidence to the contrary. 



III. ra-SITE APPROXIMATIONS 



One of the most common theoretical approaches to Markov processes with spatial struc- 
ture is a truncation of the master equation known as an n-site approximation 37|]. Such 
approximations have been applied to the TCM in efforts to determine the order of the 
transition; in this section we review and extend these results. 

The simplest method in this family is dynamic mean-field theory or the one-site approx- 
imation, in which the probability of an m-site configuration is factored into a product of m 
single-site probabilities, so that, for example, P(» • • o) ~ p 3 (l — p) where • (o) denotes 
an occupied (vacant) site and p is the fraction of occupied sites. The resulting equation for 
dp/dt yields rather poor predictions for the TCM; better results are obtained using larger 
clusters. In the n-site approximation, the equations that govern the probability distribution 
for clusters of n sites are truncated by expressing the probabilities of n + 1 site (or larger) 
clusters in terms of the n-site distribution. In the 3-site approximation, for example, we 
write P(« • • o) ~ P(« • »)P(« • o)/P(« •). As n grows, the number and complexity of the 
equations increases rapidly, but it is possible to generate the equations, and integrate them 



numerically, via a computational algorithm 

in Kef. Q, the 4-site app—o. for the TCM was found to predict a continuous 
phase transition for diffusion rates D < D t and discontinuous one for D > D t . The predicted 
value for the tricritical diffusion rate D t , however, is much smaller than that reported in 
simulations (D t ~ 0.95). Since the phase diagram predicted by the n-site approximation 
generally converges to the correct one as n — > oo, it is of interest to study the results for 
larger n. In certain cases, predictions based on a sequence of n-site approximations behave 
in a consistent manner, and can be extrapolated to provide estimates of the transition point 
and critical exponents, via the coherent anomaly method 38l. l39|. 
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For small n, the position of D t varies considerably. For example, the n < 3 approximations 
yield a discontinuous transition even for D = 0, but for n > 4 there is a tricritical point at 
some D t > 0. The estimates for D t increase gradually with n; for n = 8, one finds D t > 0.5, 
for example. On the other hand, for a fixed, large diffusion rate, the transition remains 
discontinuous, with a large jump in the order parameter, which does not diminish appreciably 
with increasing n. For D = 0.98, for example, Fig. [1] shows that all the approximations 
studied (n < 17) yield a discontinuous transition. 

Recently, Ferreira and Fontanari published results casting doubt on the utility of n-site 
approximations for the TCM 40J]. They show, for example, that for D = 0, the values of 
A c veer away from the simulation value as n is increased from 11 to 18; our studies confirm 
this observation. (We note that a nonmonotonic approach to the critical point is observed 
in a stochastic cellular automaton in which at least three particles are required for particle 
generation or survival Moreover, the values of the tricritical reproduction rate Xt(n) 
(for n < 14) appear to converge to an unphysical (negative) value as n — > oo [401 ]. 
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FIG. 1: (Color online) Order parameter versus creation rate in the one-dimensional TCM with 
D = 0.98, in the n-site approximation with n = 8, 17; points: simulation. 



Thus if the n-site approximations converge to the correct values, they do so in a non- 
monotonic fashion, such that for the cluster sizes accessible with present technology (n < 20 
or so), quantitative results for A c cannot be obtained for the TCM. Although our n-site 
approximations show stable nonvanishing gap sizes at D = 0.98 for n < 18 (Fig. [T]), we find 
that for a given level (n = 13, say), X C {D) is not a monotonic function of D, and the location 
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of the tricritical point is rather uncertain. Fig. [2] shows that on increasing the diffusion rate 
from D = 0.6 to D — 0.66, the critical point shifts to higher values, but for D > 0.66 this 
tendency reverses. On the other hand, there is no evidence of a discontinuous transition 
'or 0.67 < D < 0.71. Our analysis provides a higher A t estimate for n — 13 than found in 
4o| ; the reason for this difference is not known, since our result for A t (n = 4) agrees with 
that reported in the latter work. We observe a similar behavior for n = 14 and n = 15; our 



4Q| 



tricritical point estimates do not fit on the extrapolation line given in 

For the cluster sizes studied, Aj(n) cannot be fit with a linear function of 1/n, so that 
the n — > oo limiting value cannot be estimated with confidence. It seems likely that Xt(n) 
exhibits an oscillatory convergence with n; if so, the question of whether lim n ^ 00 D t (n) < 1, 
(that is, the existence of a discontinuous transition), cannot be resolved using the available 
n-site approximation results. 




FIG. 2: (Color online) Order parameter versus creation rate in the one-dimensional TCM in the 
13-site approximation. 



IV. SIMULATION RESULTS 

We study the TCM via Monte Carlo simulation, using several approaches that comple- 
ment earlier analyses: dependence on the initial value of the order parameter, interface 
dynamics, scaling in the presence of a weak source of activity, and scaling of the quasista- 
tionary order parameter. While some of the results would seem to provide good evidence of 
a discontinuous transition, we shall defer our conclusion until the following Section. 
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A. Initial density dependence 



If a phase transition is discontinuous, the evolution of the system should depend strongly 
on the initial condition, while at a continuous transition the evolution is toward the same 



QS state, regardless of the initial condition. In 41] . simulations of the TCM at diffusion 
rate D = 0.98 are reported, showing that at the transition (A c ~ 9.60), the value of the 
order parameter at long times depends on its initial value. For initial particle densities p(0) 
between 0.3 and unity, the system evolves to the active state, while for p(0) < 0.3 it rapidly 
approaches the absorbing state. (In studies using p(0) < 1, the initially occupied sites are 
chosen at random, uniformly over the lattice.) These results demonstrate that an active 
phase, characterized by a high value of the order parameter, is accessible starting from a 
high density, but not from a low one. The findings for high diffusion rate are in sharp 
contrast to those found for D = 0, for which the critical reproduction rate is A c = 12.015. In 
this case, the particle density attains the same QS value, starting from very different initial 
values. 

Here we extend the simulations of 41[ to much larger systems. We follow the evolution of 



the order parameter p(t) in systems of L = 2 x 10 5 sites (with periodic boundary conditions) 
for times of up to 10 9 MCS, averaging over 5 -20 realizations. For D = 0.98, we find that 
for low initial values (p(0) < 0.25), the density falls exponentially for A < 9.616 (see Figj3]). 




FIG. 3: (Color online) TCM: p(t) for D = 0.98, using various initial particle densities, and creation 
rates near the transition value, as indicated. 



For the same reproduction rates, using p(0) > 0.25, the order parameter p{t) exhibits a 
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long plateau (10 3 < t < 10 5 MCS) at a high density. At longer times p(t) decays, as expected 
in a finite system. For D = 0.5, by contrast, the order parameter curves p(t), starting from 
high and low initial values, attain a common value at long times (see FigH]). 

10' f 1 1 1 j 
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t 



FIG. 4: (Color online) p(t) as in Fig. 0J but for D = 0.5 



B. Fluctuating boundary studies 

In these studies the initial configuration consists of two blocks, one fully occupied, of no 
sites, and the other, of L — n sites, completely empty. For D = 0.98 and A = 9.60, one 
finds that in the initially occupied region, the particle density quickly relaxes to its QS value 
of about pqs — 0.83. The ensuing evolution is characterized by the drift of the boundaries 
between active and empty regions. A given realization stops either when it attains the 
absorbing state, or when the number of particles indicates that the active phase has filled 
the entire system (we use a particle number of N stop = 0.84L as the criterion for this event). 
Figure [5] shows a typical history for L = 2000 and n = 1500. (The graph shows the mean 
density in blocks of 50 sites, with time increasing downward, in steps of 10 000 time units 
between each density profile.) The boundaries between active and inactive regions appear 
to follow independent, unbiased random walks. 

If the boundaries can be represented by independent, unbiased random walkers, then the 
number N(t) of particles (which is approximated by pqs times the size of the active region), 
should also follow an unbiased random walk. The walk starts at N = uq and is subject to 



absorbing frontiers at N = and N = N stop . Well known results on random walks 42j then 
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FIG. 5: Space-time evolution of particle density for D = 0.98, A = 9.60, and L = 2000, starting 
from a fully occupied region of 1500 sites and the remainder empty. Time increases downward, 
with each sweep (at intervals of 10 000 time units) showing the density profile averaged over blocks 
of 50 sites. 

imply that the probability p stop of reaching N stop before iV = is given by n /N stop . We 
estimate p s top m se ts of 100 realizations, on rings of 500, 1000, and 2000 sites; the linear trend 
evident in Fig. [H] supports the fluctuating boundary interpretation. This in turn suggests 
that at high diffusion rates, the dynamics of active and inactive domains is effectively that 
of CDP. 

Away from the phase transition, we expect the interface to drift on the average, advancing 
into the inactive region for A > A c and vice- versa. We determined the mean interface velocity 
for D = 0.98 on rings of 3000 sites, using no = 1500. After allowing the system to relax for 
50 000 MCS, we record the density profile pi, and determine the interface position X, via the 



12 



1.0 



0.8 



0.6 



Q. 

O 

Q. 



0.4 



0.2 



0.0 



/ + 



+ 
□ 



/ □ 



0.0 0.2 0.4 0.6 0.8 



1.0 



FIG. 6: (Color online) Probability p s top versus po 
initial configuration. +: L = 500; open squares: L 
Pstop = no /N s top as expected for a random walk. 



= uq/L in the TCM with an inhomogeneous 
= 1000; filled squares: L = 2000. Dashed line: 



criterion p{xj) = Pb/2, with ps the bulk particle density at the A value of interest. We then 
allow the system to evolve for an additional 20 000 MCS, and again determine the interface 
position. Fig. [7] shows the interface drift velocity v, as determined in samples of 3 x 10 4 
realizations, varying linearly with A — A c ; linear regression yields v = for A = 9.60(1), 
in agreement with other estimates of the transition point. (We verified that the interface 
velocity obtained using a ring of 5000 sites, and observation times of 10 5 and 2 x 10 5 MCS, 
agrees to within 2% with the value obtained using the smaller system.) 

The above results show that for large D the system is divided into well defined active 
and empty regions. A typical configuration of a large system (2 x 10 4 sites) at the transition 
(D = 0.98, A = 9.608) bears this out. In order to visualize the configuration of a large 
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FIG. 7: (Color online) Interface drift velocity v versus creation rate A, for D = 0.98, L = 3000. 
The straight line is a least-squares linear fit to the data. 



system, we plot the cumulative particle number S(x) = J2f=i &i versus position x, where <jj 
is an indicator variable taking values of and 1 at empty and occupied sites, respectively. 
Thus the local density p(x) corresponds to the slope of the graph at x, with empty regions 
corresponding to horizontal lines. In these studies we initially occupy half the sites, ran- 
domly, so that initially the local density is ~ 0.5. In the initial phase of the evolution, the 
global density rapidly grows to about 0.8; thereafter it begins to fluctuate, as empty regions 
form. The configuration shown in Fig. [HI for a time of about 2.6 x 10 8 MCS (comparable to 
the simulation times in [24|), consists of a series active regions, with density p ~ 0.845, and 
empty regions, giving an overall density of 0.25. Of note is the high density in the active 
regions, and the similarity of the density in active regions separated by large inactive gaps. 
We verified that for this choice of D and A, the density in active regions p a = 0.840(5), in- 
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dependent of system size and of overall density. The system reaches the absorbing state via 
fluctuations of the boundaries between active and inactive regions, which eventually drive 
the active fraction to zero, while the active region density remains constant. For comparison, 
in Fig. [S]we also show a typical configuration for D = and A = A c (0) = 12.015. The initial 
condition is the same as for D = 0.98, and the simulation is again halted when p falls to 
0.25, which occurs at t ~ 17000 MCS for these parameters. In this case the empty regions 
are typically much smaller than under rapid diffusion, and the local density in active regions 
varies considerably. 




10000 20000 10000 20000 



X X 

FIG. 8: Cumulative particle number S(x) versus position x in a typical configuration. Left: 
D = 0.98, A = 9.608; right: D = 0, A = 12.015. 



C. Effect of a weak source 

We turn next to studies of the stationary order parameter in the presence of a weak 
source, h, defined as the probability, per unit time and per vacant site, to insert a particle. 
In [3] a weak source was used to demonstrate hysteresis; here, we consider the scaling of 
the order parameter as h — > at the transition. For D = 0, we determine p at A = 12.015, 
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varying h between 10~ 9 and 10~ 5 . For each h value, a series of lattice sizes (from 5000 to 
50 000 sites) are used to estimate the limiting infinite-size value of p. We verify the scaling 
law p oc h l ' s \ with l/6 h = 0.109(1), in agreement with the value expected for directed 
percolation in one spatial dimension, l/5h = 0.10825(3) 1]. 

For D = 0.98 and A = 9.60, we observe a very different scenario. For a given value of 
h, two values of p are found, depending on the initial density. For large initial densities, 
p approaches a value of about 0.815 as h — > 0, while for a low initial density, p oc h 
(see Fig. [9]). The order parameter, moreover, is essentially independent of system size for 
L > 1000. These results are consistent with a discontinuous transition, and the absence of 
critical scaling, for D = 0.98. 
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FIG. 9: (Order parameter versus source intensity h for D = 0.98, A = 9.60, L = 1000. Upper 
curve: initial density unity; lower: initial density zero. 



D. Quasistationary order parameter 

To close this section we report results on the QS value of the order parameter as a 
function of system size. In Ref. estimates for the limiting (L — > oo) value of pqs 
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were found to exhibit a discontinuity at the transition, for D = 0.95; here we focus on 
D = 0.98. We determined p in QS simulations of duration t max = 10 7 MCS (for L = 100) 
up to t max = 2 x 10 9 MCS (for L = 10 4 ), allowing the first 10% of the time for relaxation. 
In Fig. [TD] we plot the QS order parameter versus 1/L, for A values near the transition. 
The curves divide into two families. One set (for A < 9.605) approaches zero as L — > oo, 
while for larger values of A, the density approaches a nonzero limiting value. In the minute 
interval 9.605 < A < 9.610 the limiting (L — ► oo) value of the order parameter p jumps 
from zero to about 0.6. The inset shows that at the transition point, A = 9.6084, p decays 
exponentially with system size. At a continuous transition one expects the density to decay 
as a power-law, pq S ~ L~ l3 ^ u± . 




FIG. 10: QS order parameter p versus reciprocal system size for D = 0.98 and (lower to upper) A 
= 9.60, 9.605, 9.61, 9.62, and 9.65. Inset: semi-log plot of p versus L for A = 9.6084. 
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V. WHY THE TRANSITION IS CONTINUOUS, AND WHY IT APPEARS TO 
BE DISCONTINUOUS 



At various points in this discussion we have drawn an analogy between the TCM at 
high diffusion rates and compact directed percolation. The picture that emerges from the 
simulations reported above, and from some of the earlier studies 
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231 ]. is that in this 



regime the system contains essentially two kinds of regions, one of high density, the other 
empty. The boundaries between these regions perform independent random walks, leading 
eventually to extinction of activity. 

The above scenario corresponds to the phase transition in CDP, and suggests that we 
construct a reduced description in which blocks of i sites in the active region (with particle 
density on the order of p a ~ 0.85) correspond to sites in state 1 of the CDP, and blocks of 
i empty sites correspond to sites in state in the CDP. The dynamics of the CDP consists 
exclusively of random walks performed by the interfaces between strings of 0s and strings of 
Is. As we vary A through its critical value in the TCM, the drift velocity of the interfaces in 
the corresponding CDP passes through zero, and the asymptotic density of Is jumps from 
zero to one. 

If the above caricature of the TCM as an effective compact directed percolation model 
were valid, a discontinuous transition would be guaranteed. There are, however, two addi- 
tional processes that must be taken into account. Evidently, gaps (strings of empty sites) 
can arise within active regions, else the TCM starting from a fully occupied lattice would 
never reach the absorbing state. A fundamental point is that the process of gap nucleation, 
while essential to the TCM dynamics, occurs at an extremely small rate. By "gap nucle- 
ation" we mean the generation of a gap large enough (of g* sites, say) that its boundaries 
fluctuate independently of one another. Gaps of size g > g* are equally likely to grow or to 
shrink, whereas smaller gaps tend to shrink, due to particles diffusing in from the adjacent 
occupied regions. 

The data of Fig. 3 permit an order of magnitude estimate of the rate of gap nucleation: 
the density begins to fall appreciably from its plateau value at t ~ 10 5 , and there are C(10 5 ) 
sites in the system, giving a rate of k ~ 10 -10 per site. In Fig. [IT] we plot the mean first- 
passage time t for the appearance of a gap of size g in a system of 10 4 sites (parameters 
D = 0.98 and A = 9.608), starting with all sites occupied. (In these studies t is estimated 
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using samples of N r = 50 to 500 realizations, with smaller N r for larger system sizes.) The 
first-passage time grows rapidly for smaller sizes and then crosses over to a slower growth 
around g ~ 30, at which point t ~ 6 x 10 4 . Identifying this crossover size as g* gives a 
nucleation rate of k ~ 10 -9 per site. 

It is natural to take the block size t in the CDP mapping as the critical gap size g*, so 
that a one-site gap in the equivalent CDP process is equally likely to grow or to shrink to 
zero. Since we map g* TCM sites to a single site in the CDP, the effective gap nucleation 
rate in the latter is then n e jf = g*n. The estimates for the nucleation rate and for g* given 
above yield « e // in the range 10~ 9 — 10~ 7 per CDP site. 

10 7 
10 6 
10 5 
10 4 
10 3 
10 2 

10° 10 1 10 2 10 3 

g 

FIG. 11: Mean first-passage time t for appearance of a gap of size g. Parameters L = 10 4 , D = 0.98, 
A = 9.608. 

As noted above, a small gap can shrink when particles diffuse in from outside. Can a large 
gap be destroyed in this manner? To answer this, consider an interface (treated as fixed, 
for the sake of this argument) between large active and empty regions: the local density 
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p(x) = p a for x < and p ~ for x > 0. Think of the edge of an active region as a particle 
source. Particles are emitted at a rate of order s — 1 — D, diffuse away at rate D, and decay 
at rate 7 = (1 — D)/(l + A). In the stationary state, a continuum diffusion analysis yields 
a local density for x > of 

r oo „-x 2 /4Dt 

p(x) = s / die-* ; (1) 



giving p(x) ~ se^/" 1 with an interface width w = JD/j. For D = 0.98 and A = 9.60 this 
gives w ~ 23; for these parameters simulation shows that near the edge of a large gap, the 
density decays to zero oc e~ x ' . Nucleation of an active region inside a gap occurs at a rate 
« (1 — D)[p(x)} 3 , making the probability of nucleating activity deep within a gap negligible. 
For x = 100 and D and A as above, for example, we find a nucleation rate of ~ 10~ 13 . The 
essential point is that p(x) decays exponentially, so that nucleation of activity is limited to 
the neighborhood of the edges, whereas the nucleation of gaps can occur anywhere inside an 
active region. 

Thus we conclude that the TCM at a high diffusion rate is equivalent to compact directed 
percolation with a very small, but nonzero rate of gap formation within clusters of Is. But 
this process is in turn equivalent (insofar as scaling properties are concerned) to the Domany- 



Kinzel cellular automaton (DKCA) [15| with pi ~ 1/2 and p 2 = 1 — K e ff- (Recall that in 
the one-dimensional DKCA, pi is the probability of a site taking state 1, given one neighbor 
in state 1, and one in state 0, at the preceding time, and that p 2 is the probability of state 
1 given that both neighbors are in state 1 at the previous step.) CDP corresponds to the 



line P2 — 1. From the work of Janssen [43] and of Liibeck 44J we know that for p 2 = 1 — e, 
the transition occurs at p c = 1/2 + 0{\J~e). 

The phase transition of the DKCA is discontinuous only for p 2 — 1; for any p 2 < 1, it 
is continuous and belongs to the DP universality class. We are led to the same conclusion 
regarding the TCM: for any D < 1, there is a small but finite rate of nucleating gaps within 
active regions, so that the effective value of p 2 is slightly less than unity. 

Given the nearness of the equivalent DKCA to the line p 2 = 1, it is not surprising that 
simulations of the TCM using lattice sizes L and simulation times t m yield an apparently 
discontinuous phase transition, for Lt m < It is only for large systems and long simula- 
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;ion times that the effects of gap nucleation become apparent, as in the studies of 20j and 



24| . (Indeed, the results shown in Fig. 3 are also compatible with DP-like decay of the order 



parameter.) But numerical studies of the stationary order parameter can be expected to 
show discontinuous behavior, this near the CDP line. For similar reasons, it is not surprising 
that n-site approximations, using clusters of fewer than twenty sites, miss the effect of gap 
nucleation; much larger clusters would be needed to capture this properly. 

It is perhaps worth recalling that CDP-like spreading behavior is observed in a surface 
modified version of directed percolation in one spatial dimension [45j. In this case, propa- 
gation of activity at the edges of the active region occurs with a different creation rate (A', 
say), than in the bulk, which has a creation rate of A^. Let A c denote the critical creation 
rate for the original problem, that is, for A' = A#. For A' < A c , the phase transition occurs 
at some bulk creation rate A# > A c , which means that the bulk has a finite activity density 
even for t — > oo. Thus the active cluster is compact, and the scaling behavior is that of 
CDP. The essential difference between this model and the TCM is that in the former case 
large gaps cannot be nucleated within the active region: since A# > A c at the transition, 
gaps tend to shrink. 



VI. DISCUSSION 



We have presented various pieces of evidence suggesting that the order parameter is a 
discontinuous function of the creation rate in the triplet creation model at high diffusion 
rates (our numerical studies focus on D = 0.98). We do not find evidence of stable coexisting 
active and inactive regions; the boundaries between these regions are observed to fluctuate, as 
asserted in Ref. (2^]. The numerical evidence in favor of a discontinuous transition includes 
hysteresis (under a weak source), the initial-density dependence of the order parameter at 
later times (with or without a particle source), the high particle density within active regions, 
the random- walk- like fluctuations of the boundaries between active and inactive regions, and 
the apparent absence of power-law scaling of the QS order parameter as a function of system 
size, at the transition point. 

Despite the numerical evidence in favor of a discontinuous transition, a mapping of the 
TCM to a effective dynamics resembling that of compact directed percolation, but with 
a very small gap nucleation rate, leads to the conclusion that the transition is continuous 
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for any D < 1. The studies reported above, suggesting a discontinuous transition, were 
performed using relatively small systems and/or limited times. For example, to see the 
effect of gap nucleation in the studies with a weak source of activity, the source strength h 
would have to be much smaller than the gap nucleation rate k. 

The failure of n-site approximations to give a clear indication of the nature of the tran- 
sition may again be attributed to the very small gap nucleation rate for D ~ 1. These 
approximations, however, are problematic even for D = 0: the values predicted for A c do 
not converge monotonically to the correct value with increasing n j^ . 

The continuous nature of the transition in the one-dimensional TCM was of course as- 
serted some time ago by Hinrichsen {^(J, and received further support in Park's simulations 



241 ] . Our argument nevertheless contributes to an intuitive understanding of this result, 



and might provide the basis for a rigorous demonstration of the continuous nature of the 
transition. 
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